import pegasus as pg
Count Matrix File¶
For this tutorial, we provide a count matrix dataset on Human Bone Marrow with 8 donors. It is stored in pegasus hdf5 format (with file extension ".h5sc").
Download tutorial data
!wget https://github.com/klarman-cell-observatory/scCloud-tutorial-data/raw/master/MantonBM_nonmix_subset.h5sc
When the downloading is finished, you can load the file using pegasus read_input function:
adata = pg.read_input("MantonBM_nonmix_subset.h5sc")
adata
The count matrix is manipulated as an Annotated Data Matrix structure (see anndata.AnnData for details; the figure below is also borrowed from this link):
It has 5 major parts:
- Raw count matrix:
adata.X, a Scipy sparse matrix, with rows the cell barcodes, columns the genes/features:
adata.X
This dataset contains $48,000$ barcodes and $33,694$ genes.
- Cell barcode attributes:
adata.obs, a Pandas data frame with barcode as the index. For now, there is only one attribute"Channel"referring to the donor from which the cell comes from:
adata.obs.head()
adata.obs['Channel'].value_counts()
From the summary above, you can see that each of the 8 donors (numbered from 1 to 8) gives $6,000$ cells.
- Gene attributes:
adata.var, also a Pandas data frame, with gene name as the index. For now, it only has one attribute"gene_ids"referring to the unique gene ID in the experiment:
adata.var.head()
- Unstructured information:
adata.uns, a Python hashed dictionary. It usually stores information not restricted to barcodes or features, but about the whole dataset, such as its genome reference:
adata.uns['genome']
- Finally, embedding attributes on cell barcodes:
adata.obsm; as well as on genes,adata.varm. We'll see it in later sections.
pg.qc_metrics(adata)
The metrics considered are:
- Number of genes: by default, keep cells with $500 \leq \text{# Genes} < 6000$;
- Number of UMIs: by default, keep cells with $100 \leq \text{# UMIs} < 600,000$;
- Percent of Mitochondrial genes: by default, keep cells with percent $< 10\%$;
- Percent of cells: by default, genes that are expressed in at least $0.05\%$ of cells are marked as
robust.
For details on customizing your own thresholds, see documentation.
Numeric summaries on filtration on cell barcodes and genes can be achieved by get_filter_stats method:
stats_samples, stats_genes = pg.get_filter_stats(adata)
The results are two pandas data frames, one for samples, the other for genes:
stats_samples
stats_genes
You can also check the QC stats via plots. First is the violin plot:
pg.violin(adata, keys = ['n_genes', 'n_counts', 'percent_mito'], by = 'passed_qc')
Cells marked with passed_qc == True are those passing QC, while the others will be filtered out.
Then is scatterplot on the QC metrics on samples:
pg.scatter(adata, 'n_genes', 'n_counts', color = 'passed_qc')
This plot indicates that number of UMIs and number of genes are positively correlated.
pg.scatter(adata, 'n_genes', 'percent_mito', color = 'passed_qc')
The plot above indicates that cells with high percent of mitochondrial genes tend to have fewer signals.
Now filter data based on QC metrics set in qc_metrics:
pg.filter_data(adata)
adata
After filtration, $34,654$ cells ($72.20\%$) and $23,460$ genes ($69.63\%$) are kept.
We can also view the cells after filtration with respect to donors:
adata.obs['Channel'].value_counts()
Normalization and Logarithmic Transformation¶
After filtration, we need to first normalize the distribution of cells w.r.t. each gene to have the same count (default is $10^5$, see documentation), and then transform into logarithmic space by $log(x + 1)$ to avoid number explosion:
pg.log_norm(adata)
For the downstream analysis, we may need to make a copy of the count matrix, in case of coming back to this step and redo the analysis:
adata_trial = adata.copy()
Highly Variable Gene Selection¶
Highly Variable Genes (HVG) are more likely to convey information discriminating different cell types and states. Thus, rather than considering all genes, people usually focus on selected HVGs for downstream analyses.
You need to set consider_batch flag to consider or not consider batch effect. At this step, set it to False:
pg.highly_variable_features(adata_trial, consider_batch = False)
By default, we select 2000 HVGs using the pegasus selection method. Alternative, you can also choose the traditional method that both Seurat and SCANPY use, by setting flavor == 'Seurat'. See documentation for details.
Below is a scatterplot of highly variable genes for this dataset. Each point stands for one gene. Orange points are selected to be highly variable genes, which account for the majority of variation of the dataset.
pg.variable_feature_plot(adata_trial)
Below lists the top 5 HVGs:
adata_trial.var.loc[adata_trial.var['highly_variable_features']].sort_values(by = 'hvf_rank').head()
Principal Component Analysis¶
To reduce the dimension of data, Principal Component Analysis (PCA) is widely used. Briefly speaking, PCA transforms the data from original dimensions into a new set of Principal Components (PC) of a much smaller size. In the transformed data, dimension is reduced, while PCs still cover a majority of the variation of data. Moreover, the new dimensions (i.e. PCs) are independent with each other.
pegasus uses the following method to perform PCA:
pg.pca(adata_trial)
By default, scc.pca uses:
- Before PCA, scale the data to standard Normal distribution $N(0, 1)$, and truncate them with max value $10$;
- Number of PCs to compute: 50;
- Apply PCA only to highly variable features.
See its documentation for customization.
To explain the meaning of PCs, let's look at the first PC (denoted as $PC_1$), which covers the most of variation:
coord_pc1 = adata_trial.uns['PCs'][:, 0]
coord_pc1
This is an array of 2000 elements, each of which is a coefficient corresponding to one HVG.
With the HVGs as the following:
adata_trial.var.loc[adata_trial.var['highly_variable_features']].index.values
$PC_1$ is computed by
\begin{equation*} PC_1 = \text{coord_pc1}[0] \cdot \text{HES4} + \text{coord_pc1}[1] \cdot \text{ISG15} + \text{coord_pc1}[2] \cdot \text{TNFRSF18} + \cdots + \text{coord_pc1}[1997] \cdot \text{S100B} + \text{coord_pc1}[1998] \cdot \text{MT-CO1} + \text{coord_pc1}[1999] \cdot \text{MT-CO3} \end{equation*}Therefore, all the 50 PCs are the linear combinations of the 2000 HVGs.
The calculated PCA count matrix is stored in the obsm field, which is the first embedding object we have
adata_trial.obsm['X_pca'].shape
For each of the $34,654$ cells, its count is now w.r.t. 50 PCs, instead of 2000 HVGs.
Nearest Neighbors¶
All the downstream analysis, including clustering and visualization, needs to construct a k-Nearest-Neighbor (kNN) graph on cells. We can build such a graph using neighbors method:
pg.neighbors(adata_trial)
It uses the default setting:
- For each cell, calculate its 100 nearest neighbors;
- Use PCA matrix for calculation;
- Use L2 distance as the metric;
- Use hnswlib search algorithm to calculate the approximated nearest neighbors in a really short time.
See its documentation for customization.
Below is the result:
print("Get {} nearest neighbors (excluding itself) for each cell.".format(adata_trial.uns['pca_knn_indices'].shape[1]))
adata_trial.uns['pca_knn_indices']
adata_trial.uns['pca_knn_distances']
Each row corresponds to one cell, listing its neighbors (not including itself) from nearest to farthest. adata_trial.uns['pca_knn_indices'] stores their indices, and adata_trial.uns['pca_knn_distances'] stores distances.
Clustering and Visualization¶
Now we are ready to cluster the data for cell type detection. pegasus provides 4 clustering algorithms to use. In this tutorial, we use the Louvain algorithm:
pg.louvain(adata_trial)
As a result, Louvain algorithm finds 19 clusters:
adata_trial.obs['louvain_labels'].value_counts()
We can check each cluster's composition regarding donors via a composition plot:
pg.composition_plot(adata_trial, by = 'louvain_labels', condition = 'Channel')
However, we can see a clear batch effect in the plot: e.g. Cluster 10 and 13 have most cells from Donor 3.
We can see it more clearly in its FIt-SNE plot (a visualization algorithm which we will talk about later):
pg.fitsne(adata_trial)
pg.embedding(adata_trial, basis = 'fitsne', keys = ['louvain_labels', 'Channel'])
Batch Correction¶
Batch effect occurs when data samples are generated in different conditions, such as date, weather, lab setting, equipment, etc. Unless informed that all the samples were generated under the similar condition, people may suspect presumable batch effects if they see a visualization graph with samples kind-of isolated from each other.
For this dataset, we need the batch correction step to reduce such a batch effect, which is observed in the plot above:
pg.highly_variable_features(adata, consider_batch = True)
pg.correct_batch(adata, features = 'highly_variable_features')
adata.uns['fmat_highly_variable_features'].shape
Batch correction result is stored at adata.uns['fmt_highly_variable_features']. By default, it uses adata.obs['Channel'] as the batch key, i.e. there are 8 batches in this dataset.
See its documnetation for customization. Moreover, if you want to use some other cell attribute(s) as the batch key, consider using scc.set_group_attribute method (see details).
Repeat Previous Steps on the Corrected Data¶
As the count matrix is changed by batch correction, we need to redo steps starting from PCA down to clustering:
pg.pca(adata)
pg.neighbors(adata)
pg.louvain(adata)
Let's check the composition plot now:
pg.composition_plot(adata, by = 'louvain_labels', condition = 'Channel')
If everything goes properly, you should be able to see that no cluster has a dominant donor cells. Also notice that Louvain algorithm on the corrected data finds 16 clusters, instead of the original 19 ones.
Also, FIt-SNE plot is different:
pg.fitsne(adata)
pg.embedding(adata, basis = 'fitsne', keys = ['louvain_labels', 'Channel'])
Similarly, if everything goes properly, you should be able to see that cells from different donors are better mixed in the right-hand-side plot.
Visualization¶
tSNE Plot¶
In previous sections, we have seen data visualization using FIt-SNE. FIt-SNE is a fast implementation on tSNE algorithm. Including FIt-SNE, pegasus provides 3 different tSNE plotting methods:
pg.fitsne: FIt-SNE plot, using fitsne package. See detailspg.tsne: tSNE plot, using Multicore-TSNE package. See detailspg.net_tsne: approximated tSNE plot with speed up using Deep Neural Networks (DNN) models. See details
UMAP Plot¶
Besides tSNE, pegasus also provides UMAP plotting methods:
pg.umap: UMAP plot, using umap-learn package. See detailspg.net_umap: Approximated UMAP plot with DNN model based speed up. See details
Below is the UMAP plot of the data using umap method:
pg.umap(adata)
pg.embedding(adata, basis = 'umap', keys = ['louvain_labels', 'Channel'])
Differential Expression Analysis¶
With the clusters ready, we can now perform Differential Expression (DE) Analysis, which is crucial for identifying the cell type of each cluster.
Now use de_analysis method to run DE analysis. We use Louvain result here.
Notice that for notebooks running on Terra, you need to explicitly set temp_folder parameter to avoid the incorrect "No space left" error:
pg.de_analysis(adata, cluster = 'louvain_labels', auc = False, t = True, fisher = False, mwu = False,
temp_folder = "/tmp")
By default, DE analysis runs
auc: Area under ROC (AUROC) and Area under Precision-Recall (AUPR).t: Welch’s t test.
Alternatively, you can also run the follow tests by setting their corresponding parameters to be True:
fisher: Fisher’s exact test.mwu: Mann-Whitney U test.
DE analysis result is stored with key de_res (by default) in varm field of data. See documentation for more details.
To load the result in a human-readable format, use markers method:
marker_dict = scc.markers(adata)
By default, markers:
- Sort genes by WAD scores in descending order;
- Use $\alpha = 0.05$ significance level on q-values for inference.
See documentation for customizing these parameters.
Let's see the up-regulated genes for Cluster 1:
marker_dict['1']['up']
Among them, TRAC worth notification. It is a critical marker for T cells.
We can also use Volcano plot to see the DE result. Below is such a plot w.r.t. Cluster 1 with log fold change being the metric:
pg.volcano(adata, cluster_ids = ['1'])
To store a specific DE analysis result to file, you can write_results_to_excel methods in pegasus:
pg.write_results_to_excel(marker_dict, "MantonBM_nonmix_subset.louvain_labels.de.xlsx")
Cell Type Annotation¶
After done with DE analysis, we can use the test result to annotate the clusters.
celltype_dict = pg.infer_cell_types(adata, markers = 'human_immune', de_test = 't')
cluster_names = pg.infer_cluster_names(celltype_dict)
infer_cell_types has 2 critical parameters to set:
markers: Either'human_immune','mouse_immune','human_brain','mouse_brain', or a user-specified marker dictionary.de_test: Decide which DE analysis test to be used for cell type inference. It can be either't','fisher', or'mwu'.
infer_cluster_names by default uses threshold = 0.5 to filter out candidate cell types of scores lower than 0.5.
See documentation for details.
Next, substitute the inferred cluster names in data:
adata.rename_categories('louvain_labels', cluster_names)
And plot data with cluster names. You can either plot names on data or not, by setting legend to "data":
pg.embedding(adata, basis = 'fitsne', keys = ['louvain_labels'])
pg.embedding(adata, basis = 'fitsne', keys = ['louvain_labels'], legend='data')
marker_genes = ['CD38', 'JCHAIN', 'FCGR3A', 'HLA-DPA1', 'CD14', 'CD79A', 'MS4A1', 'CD34', 'TRAC', 'CD3D', 'CD8A',
'CD8B', 'GYPA', 'NKG7', 'CD4', 'SELL', 'CCR7']
pg.dotplot(adata, keys = marker_genes, by = 'louvain_labels')
Heatmap¶
pg.heatmap(adata, keys = marker_genes, by = 'louvain_labels')
Violin Plot¶
pg.violin(adata, keys = ['TRAC'], by = 'louvain_labels', width = 900, height = 450)
Feature Plot¶
pg.embedding(adata, basis = 'fitsne', keys = ['TRAC', 'CD79A', 'CD14', 'CD34'])
Cell Development Trajectory and Diffusion Map¶
Alternative, pegasus provides cell development trajectory plots using Force-directed Layout (FLE) algorithm:
pg.fle: FLE plot, using Force-Atlas 2 algorithm in forceatlas2-python package. See detailspg.net_fle: Approximated FLE plot with DNN model based speed up. See details
Moreover, calculation of FLE plots is on Diffusion Map of the data, rather than directly on data points, in order to achieve a better efficiency. Thus, we need to first compute the diffusion map structure:
pg.diffmap(adata)
By default, diffmap method uses:
- Number of Diffusion Components = $100$.
- Compute diffusion map from PCA matrix.
adata.obsm['X_diffmap'].shape
Now we are ready to calculate the cell developing trajectory. Below is FLE plot using net_fle method:
pg.net_fle(adata)
pg.embedding(adata, basis = 'net_fle', keys = ['louvain_labels'], width = 700, height = 700)
(Advanced) More on Clustering¶
In previous sections, we use Louvain algorithm for clustering. Including Louvain, pegasus provides 4 algorithms:
pg.louvain: Louvain algorithm, using louvain-igraph package.pg.leiden: Leiden algorithm, using leidenalg package.pg.spectral_louvain: Spectral Louvain algorithm, which requires Diffusion Map.pg.spectral_leiden: Spectral Leiden algorithm, which requires Diffusion Map.
We have used louvain method. Its documentation is here. By default, it uses resolution 1.3, does clustering via PCA matrix, and stores result with key louvain_labels in obs field of data.
Leiden Algorithm. leiden method (See details), by default, uses resolution 1.3, does clustering via PCA matrix, runs the algorithm until reaching its optimal clustering (n_iter parameter), and stores result with key leiden_labels in obs field of data.
However, with the default n_iter parameter, it can take a very long time to converge, while the result doesn't improved much. So instead, you may think of feed it a small positive number to put a hard specification on how many iterations to run. Here, we use 2:
pg.leiden(adata, n_iter = 2)
We can compare its result with Louvain's via UMAP:
pg.embedding(adata, basis = 'umap', keys = ['louvain_labels', 'leiden_labels'])
Spectral Louvain Algorithm spectral_louvain method (See details), by default, uses resolution 1.3, runs KMeans on Diffusion Map, does clustering on PCA matrix, and stores result with key spectral_louvain_labels in obs field of data.
pg.spectral_louvain(adata)
Spectral Leiden Algorithm spectral_leiden method (See details), by default, uses resolution 1.3, runs KMeans on Diffusion Map, does clustering on PCA matrix, and stores result with key spectral_leiden_labels in obs field of data.
pg.spectral_leiden(adata)
Below are the UMAP plots of the clustering results of these two algorithms:
pg.embedding(adata, basis = 'umap', keys = ['spectral_louvain_labels', 'spectral_leiden_labels'])
Save Results
pg.write_output(adata, 'tutorial_results')